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I. INTRODUCTION 



The phase structure of QCD at non-zero temperature and baryon density is important 
for the physics of neutron stars and heavy ion colhsions. The approximate phase structure 
for QCD with different numbers of quark flavors has been mapped out in various mean field 
and perturbative approximations |^ The phase structure for two masless quark 
flavors (up and down) already reveals a rich sturcture. In addition to the well known chiral 
symmetry broken and restored phases, recent investigations have revealed the possibility 
of a color superconducting phase at low temperatures and relatively high densities. The 
transition to the superconducting phase as we increase fi at zero temperature is first order. 
On the other hand, in the chiral condensation regime at zero chemical potential, the phase 
transition as we increase the temperature to the unbroken mode is second order. This 
suggests that there is a regime at intermediate chemical potentials where the chiral phase 
transition is first order. Along the line separating the broken and unbroken chiral phases 
there is a tricritical point. 

One of the most pressing experimental questions is to what extent experiments at the 
relativistic heavy ion collider RHIC can explore this rich phase structure and what would 
be the experimental consequences of having a quark-gluon plasma rather than a hadronic 
plasma following a collision of heavy ions. Since the production and evolution of the quark- 
gluon plasma in a heavy ion collision might be a nonequilibrium process, one needs to 
understand the evolution of an expanding, possibly out of equilibrium, plasma. We have 
considered a toy model, which has several properties in common with two flavor massless 
QCD to explore these nonequilibrium evolutions. The model we have found is a 1 + 1 
dimensional model of self interacting fermions, that has, in the leading order in large-N 
(LOLN) approximation a phase diagram with properties similar to that of massless two 
flavor QCD such as a tricritical point as well as a superconducting phase transition as one 
increases fi at low temperature. Since an ultrarelativistic collision leads to an essentially 
one dimensional expansion at early times, it is hoped that the rate of expansion in our 
toy model will be similar to that found in QCD so that the rate the system undergoes 
the phase transition will be similar to what would be found in a a more realistic 3 + 1 
dimensional expansion. Since this model has asymptotic freedom, the coupling constants 
run logarithmically in LOLN which is a feature shared with QCD. In this paper we will 
confine ourselves for studying the dynamics near the tricritical point in our toy model which 
for that case reduces to the Gross-Neveu model. 

One of the questions important for RHIC is whether there are unambiguous experimental 
signatures resulting from a change in the nature of the phase transition as a function of the 
chemical potential (baryon density). In our toy model, the change is the difference from a 
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first order to second order transition. In actuality this change might be the change from 
a first order transition to a crossover phenomena. Our approach is to directly determine 
the time evolution of the plasma starting from an initial point on the phase diagram above 
the chiral phase transition and watch the evolution through the phase transition. From the 
time evolving fermion mode functions one can calculate many physical quantities such as the 
current current correlation function which determines the dilepton rate as well as various 
particle correlation functions such as that for the pion. In this paper we determine single 
particle distributions functions as well as the composite particle pion correlation function to 
try to find the difference in experimentally measurable quantities when a plasma evolutions 
traversing say a first rather than a second order phase transition. For the purpose of studying 
the chiral phase transition, we can restrict ourselves to just a sector of our toy model in 
which it reduces to the well known Gross-Neveu (GN) model This simpler model allows 
us to study evolutions on both sides of a tricritical point. The exact phase structure of 
the Gross-Neveu model in dimensions 2 < 4 at finite temperature and chemical potential 
has been the subject of several recent investigations |0. Using both renormalization group 
methods, dimensional reduction methods as well as strong coupling expansions, it is thought 
that the line of chiral phase transitions in all these dimensions is either second order or 
weakly first order, which is the same situation as pertains in the leading order in large-N 
calculation. What is missing in the leading order large-N is real scatterings that could lead 
to rethermalization. Therefore, the findings of our simulations that the distributions of 
fermions and antifermions goes far out of equilibrium following the transition, might easily 
be modified by a more realistic simulation. The calculations presented here must be thought 
of as presenting the first field theory simulations at finite chemical potential of an evolution 
through a chiral phase transition with a realistic expansion rate for the plasma appropriate 
to a heavy ion collision. Future studies will remedy some of the shortcomings of this toy 
model, in that inhomogeneous plasmas will be studied as well as resummation methods 
will be used in future simulations which will still be based, however on Gross-Neveu like 
models albeit in 3 + 1 dimensions as well as in 1 + 1 dimensions. The approach we take here 
is to directly solve the evolution equations of a quantum plasma in leading order in large-N. 
A complementary approach is to study critical slowing near the critical point using ideas 
from universality and dynamical critical phenomena §]. Our interest is more in having a 
complete space time picture of an evolving quark plasma and our hope is that once we resum 
the LOLN approximation using Dyson equations and consider the 3 + 1 dimensional version 
of this model that we will be able to address issues of critical slowing down. The calculations 
presented here in the toy model already exhibit the effects of how having small rather than 
zero quark masses at high temperatures change the time period of the transition. They 
also demonstrate how one can calculate all spatial correlation functions as well as the time 



4 



evolution of the Temperature and chemical potential, and how a hydrodynamic approach 
can be quite accurate before the phase transiton. 

Following a relativistic heavy ion collision, the ensuing plasma expands and cools 
traversing the chiral phase transition. In hydrodynamic simulations of these collisions 

|lO|, |ll], 0], as well as in parton cascade models and other event genarator approaches 



O, 113, IT^, |T6|, O, IT8|, one finds that it is a reasonable approximation to treat the initial 



phase of the expansion as a 1+1 dimensional boost invariant expansion along the beam (z) 
axis. In this approximation, the fluid velocity scales as z/t. In terms of the variables fluid 
rapidity 77 = ^ In (fzf) and fluid proper time r = (t^ — z^)^/^, physical quantities such as cr, e 
become independent of t], as discussed in refs. |l^, |Tl|, 0]. Such an approach was used in our 
field theory calculations of the production and evolution of disoriented chiral condensates 
in the 0(4) a model in ref.||l9|. This approximation is valid for particles produced in the 
central rapidity region. To study more peripheral collisions a full inhomogenous calculation 
must be performed. This latter study has just started and will be the subject of a future 
paper. 

These kinematical considerations translate into the expansion being homogeneous in the 
fluid rapdity rj, which allows us to convert what would be a set of partial differential equa- 
tions for the mode functions to a much simpler set of ordinary differential equations in the 
parameter r. The LOLN approximation we will use in obtaining the field equations has 



been been discussed earlier by ourselves and others in ||T^, ^ ^ |2^, |2^, ^ and 
applied to the problem of disoriented chiral condensates in |2^. Extending the boost 
invariant simulation to 3 + 1 dimensions so that transverse distributions can be studies is 
relatively simple. 

In solving the time evolution equations for the quantum fields, the initial conditions for 
the fields are specified at r = tq, that is, on a hyperbola of constant proper time. The r 
evolving energy density and pressure are obtained from the expectation value of the energy 
momentum tensor. To discuss the production of particles we introduce the concept of 
an adiabatic number operator which is an adiabatic invariant of the LOLN Hamiltonian. 
Although our equations will be valid for arbitrary initial conditions, to study the regime 
around the tricritical point we will assume that at some initial proper time the system 
can be described by a Fermi-Dirac Distribution with given /xq and Tq in the comoving 
frame. In our simulations we will also choose the initial conditions on our mode functions to 
agree with the lowest order WKB approximation result. By choosing this initial condition 
on the mode functions, the adiabatic number operator then gives a smooth interpolation 
between the initial Fermi-Dirac distribrutions described by /iq, Tq and the final outstate 
number operators. The rest of the paper is organized as follows. In section II we review 
the equlibrium properties of the GN model at finite /z and T in the LOLN approximation. 
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Particular attention is paid to the phase diagram. In section III we derive the action in curved 
coordinates in order to discuss the evolution in terms of the parameters r and 1]. Section 
IV is concerned with renormalization and obtaining exphcitly finite evolution equations. In 
section V we discuss our choice of initial conditions. In section VI we derive an expression 
for the expectation value of the energy momentum tensor and obtain expressions for the 
renormalized energy density and pressure in terms of the mode functions of the fermion field. 
In section VII we introduce the adiabatic number operator and obtain simple expressions 
for both the fermion and antifermion interpolating number operators in terms of the modes. 
In section VIII we discuss our numerical results for the proper time evolution of the effective 
fermion mass, fi, T as well as the interpolating number densities. In section IX we determine 
the time evolution of the pion correlation function. Some of the results we obtain in this 
article were summarized and presented at a Riken workshop . 



II. EFFECTIVE POTENTIAL AND PHASE STRUCTURE 

The Lagrangian for the Gross-Neveu Model is 

c = -t^a^d,^^. - ^g' {^^'^') \ (1) 

which is invariant under the discrete chiral group: \E'j — 'j^'^i- In leading order in large N 
the effective action is 

■ +trln5-^[a] , (2) 



Seff = I d X 



where S~^{x,y)[a] = (7^(9^ + a) 6{x — y). 

The phase structure of the GN model at finite temperature and chemical potential in this 
approximation has been known for a long time |^ and is displayed in Fig. |I|. 



This figure summarizes several facts. In the GN model there is spontaneous symme- 
try breaking of the discrete chiral symmetry at zero chemical potential and temperature. 
The value of the vacuum expectation value of cr at the minimum of the effective potential 
determines and is equal to the mass of the fermion in this approximation. At zero tem- 
perature, the symmetry is restored at a critical value of the chemical potential /^c = This 
phase transition is first order. At zero chemical potential the system undergoes a second 
order phase transition to the unbroken symmetry phase as the temperature is increased. As 
a result of these two facts, at some point in the phase diagram there is a tricritical point 
which can be determined numerically to be at — = .608, — = .318. These facts can be 
ascertained by studying the effective potential in LONL which is given by 

V.fA-\T,^^) = ^[In ^-1] 
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FIG. 1: Phase Structure at finite Temperature and Chemical Potential fi. Dotted lines correspond 
to approximate analytic results described in the text. 



-ir^ [In (1 + e-^^^-'^)) + ln(l + e-^(^+^))] . (3) 

P Jo ZTT 

The integrals can be evaluated in the high temperature (and small ^) regime. Keeping the 
leading terms in the expansion one obtains 

v;„(<.^^.M) = ^(lug + ^M(^= + ^)), (4) 

which leads to the relationship: 



c - ^ exp[7 J 



(5) 



for the phase transition temperature in the regime where the transition is second order and 
^ << 1. At small ji"^ one has approximately 



(6) 



In the low temperature regime for the case a < /i we can make an approximation to the 
Fermi-Dirac distribution function that again allows us to perform all the integrals analyti- 
cally and determine an approximate analytic form for the Effective Potential. We write the 
derivative of the potential in the form: 

dV (J ^ a ("^ dk, 



— = — In — + -/ — 2-tanh 
da 277 mj n Jo E 



2T 



2T 



(7) 
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where E = \/W+~o^, and then replace the function tanh(-^^) with the straight line inter- 
polation: 

tanh(x) {1 if X > 2; - 1 if x < 2and x if |x| < 2} . (8) 

Using this approximation the integrals can be performed and V determined. The results are 
shown as the dotted curve in Fig. 1. The analytic expression is given in P|. 

When T = 0, the effect of the chemical potential is the most dramatic. In that limit 
tanh(x) = e(x), and we obtain the exact result 

dV a a fV^^ dk^^ 

— = — ln^ + -/ —Q(a-cr) 

da 27T Trio E ' 



This can be integrated to give the result that for a < fi the effective potential is given by: 



1 r 2/0 1 r/i+V/^^^^i 



Veff = -W\2 H^^^^^^ ] - 1) - - a2 + df,)} , (10) 

wheras, for a > /x the effective potential is equal to its ;U = value, namely 

Veff = ^[ln (11) 
The integration constant can be fixed by choosing VJ=//(cr = 0) = 0, which yields 

Cifi) = ^ . (12) 

At T = in the broken symmetry phase the effective mass is independent of fi and is given 
hj rrif, its value when /x = 0, T = 0. When fi"^ > m'j/2 then the true minimum is at cr = 0. 
The transition at /x^ = m'j/2 is a first order transition as can be determined by eq. @ and 
eq. (p!oD . In the toy model |^ with two coupling constants, which also has a superconducting 
phase, the first order transition takes place at the point: 

2 

/.^ = ^(l-e--^), (13) 

where 5 is the difference of the inverse of the two coupling constants of the model , namely 
(5 = ^ — When the second coupling constant k — > O"*", the toy model reduces to the GN 
model. 

In our straight line interpolation of the tanh function we obtain for the tricritical point 
which seperates the regime between the first and second order phase transitions: ^ = 
.661, ^ = .31 as opposed to the exact result 

= .608, ^ = .318. (14) 

nif rrif 
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FIG. 2: Evolution olVgff as a function of T. This is for a first order transition. 

In Fig. |l| we plot the exact numerical result for the phase diagram along with these two 
approximate results. In Figs. ^ we show the evolution of the effective potentials for the 
first order and second order phase transitions when we keep fi fixed on two sides of the 
critical value and we decrease the temperature. 



III. GROSS-NEVEU MODEL IN CURVILINEAR COORDINATES 

In order to make best use of the kinematic constraint that we are in a scaling regime where 
the fluid velocity is v = z/t, we make a coordinate transformation to the light-cone variables 
T and rj, which are the fluid proper time and rapidity respectively. These coordinates are 
defined in terms of the ordinary lab-frame Minkowski time t and coordinate along the beam 
direction z by 

z = Tsmh.ri , t = Tcosh.ri. (15) 

We shall use the metric convention ( — h) which is commonly used in the curved-space 
literature. In what follows, we use Greek indices for the curvilinear coordinates r and t], 
and Latin indices for the Minkowski coordinates z and t. To obtain the Fermion evolution 
equations in the new coordinate system it is simplest to use a coordinate covariant action 



9 



V for mu=0. 1,1=0.8 



V for mu=0.1,T=0.55 




FIG. 3: Evolution of 14// as a function of T. This is for a second order transition. 



such as that used in field theory in curved spaces, even though here the curvature is zero. 
The Minkowski hne element in these coordinates has the form 



Hence the metric tensor is given by 



g^^ = diag(-l,r2) 



(16) 



(17) 



with its inverse determined from g^^Qvp = S^. This metric is a special case of the Kasner 
metric . 

The vierbein V|f transforms the curvilinear coordinates to Minkowski coordinates, 



9,. = V^VXb: 



(18) 



where rjab = diag{— 1, 1} is the flat Minkowski metric. A convenient choice of the vierbein 
for the metric (|l^ for our problem is 



so that 



V;'^ = diag{l,r} 



K'^ = diag^,- 



(19) 
(20) 
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The determinant of the metric tensor is given by 

det V = -x/— ^ = T ■ 



(21) 



The action for the Gross- Neveu Model in general curvilinear coordinates (see pof) with 
metric g^^ is 

,2 



S = j Sx (det V) 
The covariant derivative is (see [ETI]) 

where the spin connection is given by 



(22) 



(23) 



(24) 



with Y^^y^ the usual Christoffel symbol. The label i = 1 . . . N corresponding to the SU(N) 
symmetry. For the metric eq. (\[7\} (see |^) one finds 



F = F = F = 



(25) 



The coordinate dependent gamma matrices are obtained from the usual Dirac gamma 
matrices 7° via 

r = rV^'ix) . (26) 
The coordinate independent Dirac matrices 7" satisfy the usual gamma matrix algebra: 

{7^7'} = 2r7'^^ (27) 

From the action eq. (|2^ ) we obtain the Heisenberg field equation for the Fermions, 

(7^V^ + a) = , (28) 

which takes the form 



1' 



1 



7 

+ —dr^ + a 



^, = 



2tJ t 

Variation of 5* with respect to cr yields the constraint equation: 



(29) 



(30) 



which defines the rescaled coupling constant A. Since we will be interested in having 
copies of the Fermion quantum field, the rescaled coupling constant is the relevant one for 
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discussing the large- A?" limit. The lowest order in large N (LOLN) approximation is obtained 
by integrating over the Fermi degrees of freedom in the generating functional for the Green's 
function and keeping the saddle point contribution in the integral over the constraint field 
a. One obtains that the gap equation in leading order is: 



a — 



.A 
"2 



(31) 



where we have assumed there are N identical \E'i = ip. In the scaling regime v = z/t, the 
order parameter, which is the effective Fermion mass, is independent of rj and is a function 
of r only. 

For the heavy ion collision problem, we want to solve these equations subject to initial 
conditions specified on the hypcrboloid r = tq. In LOLN, specifying the initial value of 
the density matrix is equivalent to specifying the initial particle-number density and anti- 
particle number density with respect to an adiabatic vacuum state [see below] . To complete 
the specification of the initial state, the mode functions for the Fermi field also need to be 
specified at r = Tq. 

The Dirac equation reduces to its Minkowski form if we do a rescaling 

1 



V'(x) 



and introduce the conformal time u via 



(32) 



mr = e 



One then obtains 



where 



j^'du + j^ + aiu) $(a;) = , 



a(u) = ar = —e^ 
m 



(33) 



Our assumption that the evolution is homogeneous in the rapidity variable 77 allows us 
to expand the Fermion field $ in terms of Fourier modes in the momentum A;^ conjugate to 
77 at fixed conformal time 



-ikr/T] 



(34) 



The (p^ then obey 



7°^ + i'j^kr^ + a{u) 



0^(«) = O. 



(35) 
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The superscript ± refers to positive- or negative-energy solutions with respect to the adia- 
batic vacuum at r = Tq as we shall show. It is convenient to square the Dirac equation by 
letting: 



where the momentum independent spinors x"^ chosen to be the eigenstates of ij' 



(36) 



and obey the normalization condition: 



XrXs ^rs 1 



(37) 



(38) 



with r, s = ±. An explicit representation for these spinors is found in the appendix. The 
quantities f^{^)x^ ^"^^ fk{'^)x^ the two linearly independent solutions corresponding 
to the positive and negative energy solutions. When A;^ = we have that (in what follows 
the T) in A;^ when fc^ is a subscript is suppressed for notational simplicity) 



d J 
du 



fk=o = , 



so that 



The mode functions obey the second order differential equation: 



where now 



(39) 

(40) 
(41) 



If we impose the canonical anti-commutation relations on the fields 

and assume the usual canonical anti-commutation relations on the Fock space mode opera- 
tors, 

{h,bl}^{d,,dl}^2nS{k-q) , (42) 
then we obtain the condition: 



(43) 
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Taking the trace yields the normahzation condition: 



= 0, Vfc = 1 . (44) 

Using eq. (pSD one obtains 

^[0lV.] = O, (45) 

where a,b = +, — , for every a, b. Reexpressing this in terms of the f"' we have for each 
Fourier mode at all u 



fl^f^ + ^Ifl^f^ ± {f^fl^ - fl^f^) = 1 • (46) 
We also have that 

= K {ft^fk - ft^fk) ■ (47) 
Since the right hand side of eq. (^Tj) is proportional to the Wronskian and is thus independent 
of the conformal time, if we initially choose 

= irk - /r/,- , (48) 

then the two solutions 0"'", 0~ remain orthogonal at all times. This relation as well as 
eq. (^0[) can be satisfied by having 

/*+ f— 

Jk — Jk ■ 

These results can be summarized by 

= , (49) 

with a, b taking on either + or — . For both renormalization purposes as well as to in- 
troduce the concept of adiabatic number operators, it is useful to also have a WKB-like 



parameterization of the positive-energy solutions as discussed in Ref. ||T2 



f^{u) = exp I r (~^Q,iu') - du'] . (50) 

The u derivative is given by 

/+(n) = (-^(]fc- A,) /+(«), (51) 

where 



The Clk obey the real equation 

^liu)-nl, (53) 



^fc + q"fc _ ((T + ^k){0' + 3^7fc; _ 
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which is the starting point for a WKB expansion of the mode functions. Qk and Qk can be 
determined from the mode functions as follows: 



Qk = — Im 



Jk 



- 2Re 



' k fk 2f2fc 

Using the normalization condition eq. (^6|) we can show that 



Zk{u) = I/; 



±|2 



-1 



(54) 



(55) 



The normalization of the wave function, A^^ is time independent and can be evaluated at 
u = 0. Using eq. ( pS] ) and eq. (|50| ) we obtain 

Nl = 2fifc(M = 0)Zkiu = 0) . (56) 

We can now obtain an expression for the the gap equation in terms of the mode functions. 
Using the mode decompositon and the definitions: 

<b'^{k)b{q)> = 2-n5{k - q)N+{q) , 

< d\k)d{q) > = 2TT6{k - q)N4q) , (57) 



we obtain 



A 

-2 — 

2 



27 



27T 



(58) 



(1 - 2iV+(A;))0+V0^ + (2iV_(A:) - 1)0, 

Later we will choose the initial state number densities N±{k) to be Fermi-Dirac distributions 
at a given /i and T. We remark here that the k in N±{k) refers to A;^ canonically conjugate 
to rj. When we calculate the expectation value of the energy momentum tensor to identify 
the energy and pressure in a comoving system we will find that the N±{krj) corresponds to 
the comoving number densities pertinent to relativistic hydrodynamics. After some algebra 
we find 



Rk{u) 



^k = I + - ki)\m' + ^^ {ftfr - fk^f'k 



\fk 



nl~Col + 2a{a + Q,) + {^^f 



(59) 



Eq. (pq) reduces to the vacuum expression for the gap when everything becomes time 
independent so that 



[(2u;fc)(u;fe + (x)] \ 
a + r^fe, 2 



(60) 



2fii 



2cr(cr + cjfc) 
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which leads to the vacuum gap equation: 



ao = xf^ ^ = X[^ ^. (61) 



In the time evolving case we obtain instead 

A /• dk 



a = - [-^{l-N4k)-N.{k))R,{u) , (62) 

T J ZTT 



which in dimensionless form can be written as 

dk. 



a =\f^{l-N^{k)-N^{k))R,{u). (63) 

J ZTT 

We can simplify eq. ([59| ) for R{u) by using eq. (|46| ) to obtain 

i?,(n) = l-2fcj |/+(«)p. (64) 



IV. RENORMALIZATION 

In this section we will show that the renormalization of the charge A in the vacuum sector 
is sufficient to render the equation for a finite. First let us remind ourselves that the effective 
potential for the Gross-Neveu Model can be written as: 

\/ = ^ = ^-^Trln(7-9-cr), (65) 

where here the trace is only over the Dirac spinor indices. The renormalized coupling 
constant defined at arbitrary a is given by 

d'V _ I _l fdp 1 ^ 1 

dada Xn{cr) A i 27r + a' n ' ^ ^ 

It is useful to define the logarithmically divergent integral 

""^'"'^ - C'ivwrw- '''' 

In particular if we choose the renormalization point to be at the minimum of the potential, 
which also defines the fermion mass in this approximation, one has 

— U=™,=0^-=/ /^== = SmJ =— In— . 68 

da ^ A J -A 27T Ip2 -I- 77j2 277 rrij 

From the above equations, one deduces that 

XR{mf ) = n , (69) 
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and also that 



Sfmi 



As derived earlier, for the evolution problem a is given by: 

dk, 



-ln(- 



f 



a 



A 



2tt 



N^k) - N.{k))Rk{u) 



(70) 



(71) 



The apparent logarithmically divergent part of the momentum integral which we call cr'^*" 
comes from the term: 



X 



[l-2^iZ,{u)] . 



(72) 



The logarithmic divergence can be isolated by doing an adiabatic expansion of the in- 
tegrand in the expression for 5"'^™. The first order adiabatic expansion of the equation for 
the generalized mode functions Vtj. is obtained from the expression for f2| in eq. ( ^31) by 
replacing fi^, Vt^ by cD^,, cj^ in the left hand side of that expression. We further use the 
expression for ojk in terms of a expressed in eq.(^) to get: 

21 



i2 

o— + — 

UJk UJk 



(73) 



At large momentum, we therefore obtain the expansion: 
1 



1 



all 



a 

oJk 



^2 



1 



^k 



UJk J 



+ 



l2 
O 



a 

^k 



1 + ?^ 

^k 



It therefore follows that 



X 





■ a ' 


1 27 


■ i^k- 



(74) 



(75) 



We can renormalize the equation for a be appropriately subtracting this quantity when we 
renormalize the coupling constant, or we can recognize that this divergence is only apparent 
once we utilize the gap equation for the vacuum sector: 



2vr Jk^^ + mj 



dk 



1 



'^T^Jk^ 



m 



(76) 



That is we just to use the mode sum version of eq.(^) in place of A in eq. (^) and keep 
enough modes in the numerator and denominator until the answer is independent of the 
number of modes. This approach has the advantage that it also allows one to verify that 
the coupling constant is flowing according to the continuum renormalization group flow as 
we increase the cutoff. 

To explicitly renormalize the gap equation, we consider the quantity: 

_] = _a{u) ln(-,^^ 



Aij(o-) X[i{mf] 



1 ir"^ r dk 

-a{u) ln(-^) = j ^[Rf k{u) - {N^{k) + N_{k))Rk{u)], (77) 
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where Rfkiu) is given by 

Rf^[u) = Rk{u) - ^ . (78) 
From the above discussion of the divergence structure, this equation is manifestly finite. 

V. INITIAL CONDITIONS 

To solve for the proper time evolution of this system one has a to solve eq. (3.26) 

as well as the eq. ( |77D for the order parameter 

-a{u) ln(— ) = / -^[Rf ,{u) - {N^{k) + N_{k))R,{u)] . 

ZTT nijr J ZTT 

To solve eq. ( |77D requires knowledge of the initial distribution of fermions and an- 
tifermions N^{k) and N_{k). In general the only conditions needed on the distributions 
N^{k) and N-{k) are that they lead to finite number and energy densities. The condition 
on the mode functions for the energy density to remain finite at all times is that the high 
frequency modes must converge to their free field vacuum values. Since the low frequency 
modes do not effect the renormalization they can be chosen arbitrarily. Here we will choose 
the initial modes to match up with those found in the WKB approximation initially. This 
will allow us to introduce an adiabatic number operator which has the added feature that 
it will interpolate from the initial distribution of N±{k) to the final outstate distribution 
functions. Other choices are perfectly acceptable but would lead to a jump in the number 
of particles from their initial values right after the initial time. Since we are interested in 
exploring the phase diagram on both sides of the tricritical regime, we will choose our initial 
state to be described by a value of N±{k) consistent with local thermal and chemical equi- 
librium described by the parameters T and /i (here the plus sign corresponds to the Fermion 
distribution function and the minus sign to the antiFermion distribution). 

where 

u;,, = E= + a\u = 0) = ^^k^ + a^{u = 0) . 

In LOLN approximation, the initial value problem is totally specified once N± and fk and 
fk are given. The dynamics of the evolution is such that the gap equation is obeyed at all 
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times. If we start in the regime where symmetry is broken, then the initial value of the 
Fermion mass obeys the equation 

0) 



0) 



-To 



where a satisfies the gap equation eq. ([77| ) with u = 0. If we start our simulation in the 
unbroken mode then we have instead 



aiu 



0) = 



We notice that in LOLN when a = the Fermions evolve as if they are free and massless. In 
the full theory, one obtains non-trivial dynamics for the unbroken phase through fluctuations 
that are ignored in the LOLN approximation. Thus in LOLN we need to think of the a = 
case as a limit of taking initial conditions with small but finite initial mass in order to have 
nontrivial dynamics. As we will see below, finding the limit of zero mass is numerically 
well defined once initial masses are < 10^^ rrif. For our simulations an initial value of 
(t(0) = lO^'^nif is used. 

It is convenient to choose the initial tq = ~: and measure the proper time in these units. 
The adiabatic initial conditions on the mode functions / which we will discuss in more detail 
later correspond to 

ClkiO) = CukiO); 4(0) = CokiO) = , (80) 



so that using eq. 



ID we have 
/fc(0) = 



-lUJk - 



A^o) 



2c:^fe(o) 



2cu2(0) + 2cDfc(0)a(0) 



a^(0)+d;fc(0) 



(81) 



We now show that a{u = 0) = is required if we choose adiabatic initial data for the 



mode functions. For adiabatic initial conditions 

a + ujk{0) 



^(0) = A / ^iV2 



2n^'''uk{0) 



UJk 



1 - N+{k) - N^{k) 



(82) 



Eq- (|T7|) and eq. ([8^ ) are two equations for the intial conditions a and a in terms of n 
and T. However we realize that the integral for o"(0) is proportional to o"(0) times the bare 
coupling A times a finite function of cr(0), 5-(0), /i, T. Since the bare coupling goes to zero 
with the cutoff, the only value of 5"(0) consistent with the adiabatic initial conditions is 



^(0) = 



(83) 
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It then follows that 



ln^^r^-^[2- tanh ^ - tanh ^] . (86) 



so that the equation for cr(0) simplifies to 

'■(0) 

This is equivalent to the gap equation arising from the Effective Potential at finite chemical 
potential and temperature 

-oo OJk 

if we choose N±{k) to be equilibrium Fermi-Dirac distributions. 

Summarizing, for our choice of adiabatic initial conditions our mode functions initially 
are: 

Nu 

/+(0) ^-iu,f^{0). 

(87) 

We are interested in studying evolutions with and without phase transitions, as well as 
comparing the effects of traversing first vs. second order phase transitions, including the 
special case of traversing the tricritical point. We have thus chosen four separate illustrative 
starting points on the phase diagram of Fig. 1 for our numerical simulations, all assuming 
initial local thermal and chemical equilibrium. For the case of no phase transition, case (1), 
we have chosen the starting point (//q = -2; To = .3). For typical initial conditions for which 
the second order phase transition is traversed (case (2)), we choose = .5;To = .5. For 
the initial conditions /Iq — .6;Tq = .32 (case (3)) the system traverses the tricritical point. 
Finally, for a typical case where the system undergoes a first order phase transition, case 
(4) we choose (//q = -8; Tq — .3). 

We will find that before the phase transition, the system can be described in terms of a 
number distribution with a proper time evolving temperature and chemical potential. How- 
ever, after undergoing a phase transition, the adiabatic single particle distribution functions 
for Fermions and antiFermions (defined below) are far from equilibrium and cannot be de- 
scribed by a chemical potential and temperature which are independent of the momentum. 

In performing our numerical simulations, we place the system in a box of dimensionless 
length L — mfL, and choose antiperiodic boundary conditions for the Fermions. That is we 
let 

^ k^n] = ^J^tz^ , (88) 
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with n — —N, ■ ■ ■ N. In our simulations, for the system to be in a regime where the couphng 
constant flows according to the renormahzation group, one needs N — 5000. The number of 
modes increases if we want to study the reaUy long time behavior of this system. However, 
at these long times we expect hard processes which allow rethermalization, and which are 
neglected in this mean field study to become very important. 



VI. ENERGY-MOMENTUM TENSOR 



To understand the hydrodynamical properties of the evolution of the plasma, we need to 
evaluate the expectation value of the energy momentum tensor in the initial density matrix. 
Because (Tj^i,) is diagonal in the 77, r coordinate system, we can read off the comoving pressure 
and energy density from its diagonal entries. The energy-momentum tensor is defined via: 

T = (89) 

Scaling out the factor of N, and using the equation of motion we have 

^ = |^7(mV.)^' - 1^1,^%)^ - 9,u^ . (90) 

Here the parenthesis means keeping both terms in the symmetrization in ni>. If we rescale 
the fields = y/rip and use the Fourier decomposition for the rescaled fields we find for the 
expectation value of the unrenormalized T^,- 



+ (1 - 27V_(/.))(0-t9„0^ - . (91) 



Expanding in terms of the mode functions we obtain: 



e r = 



N 

^2 cA dt 



k'lo ^(l-^+(^)-^-(A:))[2a + 4Q,(a;^-cr^)|/+n, (92) 



where 



The energy density contains an infinite (quadratically divergent) cosmological term: 

'■^ dk. 



9ooK = - [ — k , 

Jo TT 



that needs to be subtracted by hand. The remaining logarithmic divergence is eliminated 
by coupling constant renormahzation. The divergence structure can be analyzed using an 
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adiabatic expansion of Qk in terms of Uk and recognizing the high momentum behavior of 
Clk is given by 

(93) 



so that 



= 7^ r[l + 0( — )1 



(94) 



Thus the divergent terms in {T^-t) are exactly the terms that appear in the effective 
potential. Namely 

^2 /rj^drv^ ~2 



N 



(95) 



Making the rescalings: 



o = GT] A = Ar; /c^ 



kr 



one recovers the unrenormalized effective potential for the Gross-Neveu model: (See eq. 
(2.21) of i). 

where here we have subtracted the (infinite) a independent cosmological constant term. 
Next using the fact that the bare coupling and the cutoff and the physical Fermion mass of 
the vacuum theory are related by the gap equation: 

dk 1 



1 
2A 



27r./A;2 + 



m 



(97) 



we obtain after subtracting the cosmological constant term 



\ sub TTI 

N 



'■A dk 



= _ 2VA;2 + a2 + 2k] = ^[In ^ - 1] + 0(-^) . (98) 

f •' 



a 



2tt /P + m 



Here we have chosen the zeropoint of the Effective Potential to be zero at the maximum 
a = 0. Thus in the true broken symmetry vacuum a = 0, we have 



V[a = mf] = - 



m 



f 
Air 



(99) 



Thus we expect (and will find) that when the system goes through the phase transition into 
the broken symmetry phase, that the energy density will approach this negative true vacuum 
value. 

The manifestly finite expression for e{T)T'^ is 



e(T]T 



27 



a' 



k'l + m 



^2 

f 



+ 2{k^-a)+Ank{a'-0i)\f, 



+ -^(N4k) + N4k)) 2a + Ank{Coi-a')\M 
Jo Ztt l 



(100) 
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We have for the expectation value of the unrenormahzed T, 



vv 



N 



0"^ i 
'2X^2 



+ (l-2iV_(fc))0-SV0fcJ ■ (101) 
Expanding in terms of the mode functions we obtain: 

Multiplying by and keeping the lowest order in the adiabatic expansion as before we 
obtain that the divergent part of the pressure is given by 



Po 



A dk 2^ 



27r v^FT^ 



A2 

'2ti 



a 
47r 



In 



cr 

4A2 



This is to be compared with divergent part of the energy density given by: 



eo = -2 



A dk 



^TFT^ = _^ + ^ ( In 
27r 27r An \ AA^ 



a 



- 1 



(103) 



(104) 



As we discussed in the momentum cutoff A acts as a noncovariant point splitting 
regulator, giving terms in the regulated {T^^) proportional to 5^5^ in which the spatial 
directions i,j = 1, 2, 3 are distinguished. Since these terms do not appear in the fi = u = 
time component, the energy density has the correct A dependence and requires no correction 
to make it agree with covariance. However, because covariance requires that the cosmological 
term is proportion to g^i,, the correct regulated pressure must satisfy 

p'o = -eo • (105) 
We need to enforce this condition by hand by adding the difference 

Po -Po = -eo - Po 

to Po which corrects for the noncovariant term induced by our momentum cutoff. We then 
need to subtract off the correct cosmological term to eliminate the quadratic divergence. 
As in the case of the energy density, the apparent logarithmic divergence gets cancelled by 
coupling constant renormalization. We then find: 

r2 ^2 ^2 A 2 

^2" 



, a" cr" (J^ ^ , A^ 

Pn T = In — T — 1 H • 

^° 2A An^ m) ^ 2n 



(106) 



Thus the renormalized and "covariantized" expression for the pressure is: 

dkr, 



I 2 

p r 



+2 



2n 



kl + d^ 



N4k) - N^k)) A{a + n^w - ^iM 
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^2 



a 



k?. + 



~ 2 
rrif 



(107) 
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A. Local Equilibrium Hydrodynamical Picture 



In this subsection we will develop a simple local equilibrium hydrodynamical model with 
an ultrarelativistic equation of state which we can compare with our field theory simulation. 
This type of model was first put forth by Landau |0. In the hydrodynamical model, Landau 
assumed that that the fiow of energy and momentum following an ultra-relativistic high 
energy collision of protons or heavy-ions behaved as an ideal fiuid fiow into the vacuum 
with initial conditions connected to a highly Lorentz contracted disc of matter. A related 
idea due to Bjorken [jlT| was to assume that because of the fiatness of rapidity distributions 



in heavy-ion collisions, the hydrodynamic fiow should posess invariance under boosts. In 
either case one has approximately that the fiuid velocity Vz = z/t and all variables only 
depend on the fiuid proper time r and are independent of the fiuid rapidity t]. In this 
hydrodynamic approach the dynamics are incorporated into the local equilbirium equation of 
state p = p(e). For a one-dimensional relativistic hydrodynamic fiow the energy momentum 
tensor is assumed to be that of an ideal fiuid 

Here a, (3 correspond to Minkowski coordinates and 

u 



V 



2 



The covariant conservation law of energy and momentum is 

pa; 

■;/9 



= . (109) 



Introducing the thermodynamic relations: 

e + p = Ts; de = Tds , 

and assuming the scaling law 

V = z/t , 

one finds for an equation of state of the form p = c^e that energy momentum conservation 



leads to the two equations: 



d^lnT = , 



dulnT + cl = 0. (110) 



Thus for an ideal one-dimension fiuid in the scaling regime one obtains: 



To V' ' 

- = (^)i+^o . (Ill) 
Co r 
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The ultrarelativistic limit has Cg = 1 (the speed of hght in our units where c = 1) for a one 
dimensional fluid. For that case 

- = (^r; ^ = (^). (112) 

We will find from our numerical simulations, that if we start in the massless (unbroken 
symmetry) regime then indeed p = e and this fall off pertains until one goes through the 
phase transition, after which the system is no longer in local equilibrium. 

When there is also chemical equilibrium, it can be shown fS^ that for a relativistic fluid 

^ = ^ , (113) 



so when that occurs we also have 



f = (^). (114) 



B. Local Equilibrium equation of state and Single Particle Distribution Functions 



We can use the equations for the renormalized energy density eq. ( p,00| ) and pressure eq. 
( p,07|) to determine the equihbrium equations for p{fi,T) and e{fi,T) and thus the equation 
of state p = p{e) when one is in chemical and thermal equilibrium . In equilibrium, a{fi,T) 
is zero if we are above the phase transition in the unbroken symmetry phase. If we are in 
the broken phase, it is instead given by the solution of the renormalized gap equation eq. 
(^). The mode functions are given by: 

IMV = _ ,\ , . (115) 

In the unbroken symmetry regime, where a = 0, this simplifies to 



\fk\-r 



2„2 _ 

2k^ 



For the unbroken symmetry case, the equations for p and e become the same {p = e) and are 
given by: 

— 2k [N4k,i^,T) + N.{k,fi,T)] , (116) 

where and are given by eq. (^9|) and we have set a = 0. In order to get a quantitative 
estimate of what will happen in our numerical simulations, we will assume that the system 
stays in local equilibrium as it evolves with T evolving as predicted by the hydrodynamical 
model so that T/Tq = tq/t. We also find from our numerical simulations that when the 
initial chemical potential is below the tricritical value the chemical potential falls similar to 
the case of chemical equilibrium: i.e. /x//io = tq/t . Let us now show that when we are 
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0.4 



0.5 




FIG. 4: From left to right: r Invariant thermal distribution A^_|_ vs. /c^ for initial conditions 
pertaining to traversing the tricritical point; r Invariant thermal distribution N- for the same 
initial conditions. 

in local chemical and thermal equilibrium N±{kri,T) becomses independent of r. Because n 
and T scale as -, and /c = — , we find that the distributions for N4- is 



exp( 



To ' 



(117) 



so if we plot these distributions against kjj they are independent of r and just depend on the 
initial values Ho and Tq. A plot of N± vs fc^ for the initial conditions of case (3), /io = -6 and 
To = .33, the initial conditions for passing through the tricritical point is shown in Figs. |. 

From this discussion, we see that a deviation from the local equilibrium hypothesis will 
show up as a change in the adiabatic distribution functions for A^±(A;^,t) from its initial 
value as the system evolves in u = logr. This deviation from will be greatest when the 
trajectory undergoes a first order phase transition. The behavior of the plasma before the 
phase transition is very well described by hydrodynamics with equation of state p = e, 
since it is basically the evolution of a non-interacting relativistic gas initially assumed in 
equilibrium and described by: 



1 r dk 



118) 



This evolution in the massless phase is shown in fig. ^. which displays the expected 
falloff. 



VII. ADIABATIC NUMBER OPERATOR 

One important question for RHIC physics is how the quark and antiquark distribution 
functions change in time since these time evolving distribution functions enters into calcu- 
lations of the particle production rates for both pions and dileptons. Of especial interest is 
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FIG. 5: From left to right: evolution of the energy density as a function of u in the massless phase; 
same evolution for the log of the energy density as a function of u. 

how these distributions for Fermions and anti-Fermions gets modified as the plasma traverses 
the chiral phase transition. Since particle number is not conserved during the evolution of 
this system, one looks for a quantity that interpolates from the initially given distribution to 
that which obtains once interactions have ceased. For this problem because of the expansion 
into vacuum, interactions are automatically diluted at late times and one eventually reaches 
the broken symmetry vacuum state. LOLN is a mean field theory, which is related to a 
field theory for a fermion with a proper time varying mass which at late times reduces to a 
free field theory. For this reason, it is possible to introduce various time dependent number 
operators which at late times reduce to the exact out state number operator. The approach 
which we follow here, namely to define adiabatic number operators, is essentially the same as 
that used when one has the problem of quantum fields in time evolving curved background 
spaces (see for example ref. |]5D[). The basic idea is to use the WKB approximation to define 
a class of adiabatic vacuums upon which to define time evolving number operators. In the 
curved space problem, the curvature automatically introduces a time evolving mass term for 
the quantum fields. In our problem the time evolving mass term must be self-consistently 
determined. 

Within the WKB approach, there are several versions of the interpolating number oper- 
ator which differ by higher terms in the WKB expansion of the equation for the generalized 
frequencies Q, which enter in a WKB parameterization for the exact mode function f^. We 
will describe the zeroth and first order number operators below. Once we determine the 
time evolving interpolating number operators, we fit them at each m to a Fermi-Dirac distri- 
bution. We will find that before the phase transition it is possible to define a temperature 
and a chemical potential which are relatively independent of momentum. However, once 
the phase transition is traversed, the interpolating number operator does not at all resemble 
an equilibrium distribution function. To define the interpolating number operator we use a 
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complete set of wave functions (m) which are related to the solution of the exact mode 
function equation at some order in the WKB expansion. The WKB expansion is detemined 
from the second order differential equation (eq. (3.39)) for the generalized frequencies ^2^- 

In general, if we introduce a new basis y^{u) which are complete and orthonormal, but 
which do not satisfy the Dirac equation, then the number operators themselves become time 
dependent and the expansion of the quantum field becomes 

Hx) = I ^[aik,u)yUu) + c^ik,u)y^iu)]e^'^^ . (119) 

This expansion is an alternative to the expansion in terms of the initial time creation and 
annihilation operators 



For this new expansion as well as the new creation and annihilation operators to obey the 
canonical anticommutation relations, the mode functions must satisfy the orthonormality 
conditions: 

yt\u)yliu) = S^' , (120) 

for a, 6 = ±. The two sets of creation and annihilation operators are related by a transfor- 
mation which preserves the canonical structure, namely the Bogoliubov transformations: 

a{k,u) = ak{u)b{k) + Pld^k) , 

c^{k,u) = -(3k{u)b{k) + ald^{k) , (121) 

with the condition 

\ak\^ + \Pk\' = I . (122) 



Inserting the Bogoliubov transformation into eq. (|119|) and identifying terms we obtain 



4>t{u) = ak{u)yt - Pk{u)y^^ , 

<p-k{u) = Pl{u)yt + aliu)y]: . (123) 

We can project out the Bogoliubov coefficients using the orthogonality of the y or the 0, 
namely: 

c^k = ytUt; PI = ytUk ■ (124) 

If we choose our initial conditions so that y = (p, then initially 

ak = 1; Pk = 0. (125) 
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For that choice, the adiabatic particle number density will agree initially with the initial time 
number density. The interpolating number operators for the Fermions and anti-Fermions 
are defined by 

N+(k,u) = {a^ {k,u)a(k,u)) , 

N_{k,u) ^ {c^{k,u)c{k,u)) , (126) 

where the expectation value is taken in the initial density matrix parametrized by and 
defined earlier. In terms of f3 we find that 

N^{k,u) = N^{k) + [1 - N+ik) - N4k)]\Pk{u)\' , (127) 

so that the total number of particles minus antiparticles is conserved. Since f3{u = 0) = for 
adiabatic initial data, and at late times we expect a to be independent of u so that N{k, u) 
becomes the out number operator at late u, the N{k^ u) interpolate between the initial and 
final values of the average phase space number density of particles. 

We also sec that if there are no particles present initially, then |/?(A;, m)^ gives the particle 
spectrum, and its derivative is related to the rate of pair production. When particles are 
present then the presence of these particle inhibits further production because of Pauli 
blocking. 

A. Zeroth Order Adiabatic Number operator 

The zeroth order in WKB wave functions are obtained from eq. (3.36) by ignoring all 
derivatives in eq. (3.39). This yields 



= ^^^^"^{^ '^^k{u')du'^ , (128) 



One easily verifies that introducing the basis functions: 

yt = Uke~'^ yZ = v^e'^ ''''^ , (129) 

v-k = I /: = X , (130) 
^2uJk{uJk + (7} 

that the spinors are orthonormal: 

u\vk = 0; u\uk = vlvk = 1 , 
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which guarantees that the orthonormahty condition eq. (7.2) is satisfied. Using the rela- 
tionship 

/9* = ytU'k , 



one finds that 



with 



2uk{Cok + &) [nl + ul + 2Qk^ + Al] 



and Afc given by eq. 



B. first order WKB interpolating number operator 



131) 



If we keep terms up to and including all first derivatives in eqs. (3.36) and (3.39), we 
obtain the first order WKB wave functions given by: 



The gk obeys the first order differential equation 



du 



where An is 



We decompose $ as follows: 



A 



fcO 



2ujk 



where now the are given by 



Vkiu) = Mu) 



d 

du 



One can then verify that the ?/^(m) are orthonormal providing we choose 



.2/ N 2u;a: 

"^kW = ^~ /- TT . , , , exp 

2uJkiuJk + (t) + AUu) 



At time m = we have A = and 

Aliu = 0) 



2uJk 



2ujk{uJk + cr) 



du — 

UJk 



(132) 



(133) 



(134) 



(135) 



(136) 



(137) 
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so that the exact and adiabatic wave functions match up. Again using the relationship 
one finds that 



with 



|«fcP = l-|/3fcp . 



and Afc given by eq. |52| and Aqa; given by eg - 11341, so that 

A A _ ^ ^ 

Afe - Afco = 7;^ ^ 77^ • 



(138) 



VIII. RESULTS OF NUMERICAL SIMULATIONS 



We have solved the simultaneous equations eq. (^Of ) and eq. (|63|) numerically by dis- 
cretizing the Fourier modes in a box of dimensionless length L using antiperiodic boundary 
conditions for the Fermion modes. Our initial conditions were described earlier and are 
based on adiabatic initial conditions and an equilibrium value for N±. We have varied the 
time step, length of the box, as well as the number of modes until the answer was insensitive 
to our choices. The sensitivity to some of these parameters will be displayed below. Since 
the phase transition occurs near m = 2 or r = 7.4, we typically continue our calculation until 
M = 4 or r = 54. In that regime of m, it was sufficient to choose L = 500 and keep 5000 
modes. The time step needed for these values was du = 0.00004. We use a fixed grid in 
the dimensionless momentum fc^ = kr of 5000 points. This was sufficient to insure that the 
range of integration in the calculation of a included physical momentum k at least of the 
order 10m f. Because of our fixed grid in /c^ = kr, we had to increase the number of modes 
we include in our evaluation of a as the proper time increases. 

In what follows we discuss some of the finite size effects of our grids. Because of the 
exponential dependence of r on m the number of modes needed for an accurate answer at 
late times grows rapidly. We can see the dependence of our answers at late proper times 
on the parameters L, the number of modes as well as the time step Am in the figures 
below. If we keep the number of modes fixed at 5000, and decrease L then we increase the 
momentum range in our integrals and improve the result for a at late times. This is seen in 
fig. § 

In fig. 1^ we show how the evolution of a depends on the time step du. We see that once 
we have a time step du = 0.00004 then our results are insensitive to any further reduction 
in time step. 



31 




FIG. 6: Evolution of o" as a function of u for two different values of the dimensionless box length. 




FIG. 7: Evolution of cr as a function of u for different values of du. 



In fig. 1^ we see how the approach to the continuum value of ^ = 1 depends on the 
number of Fourier modes. We see that at late times there is still some depedence of the 
asymptotic value on the number of modes. 

Another issue we would like to address here is the dependence of the evolution on the 
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FIG. 8: Evolution of a as a function of u for increasing numbe of Fourier modes. 



initial conditions chosen. We have seen that in LOLN, if we start with the fermion mass 
exactly zero, then the theory is non-interacting. Thus we must consider the massless theory 
as the hmit of the massive theory, with the understanding that in higher order the fluctua- 
tions ignored in LOLN will make the equation for the fermion modes nontrivial even in the 
unbroken symmetry phase. To see that the theory actually approaches a limiting behav- 
ior, we considered initial condition where either the mass is small and not zero or the time 
derivative of the mass is non zero. The limiting theory is seen to be more readily accessed 
by choosing & = and gradually letting a ^ 0. For — < 10^^ the transition point changes 

nif 

little. We illustrate this graphically below. Figure ^ displays the time evolution for different 
initial values of &. Figure [1^ displays the time evolution for different initial values of a. As 
we go to bigger initial values of cr or a the transition is more a gradual crossover rather than 
a sharp transition. 

Based on this discussion, in the simulations presented below, we will canonically use the 
values du = 0.00004, L = 500, N = 5000, and a = 10~^. 

A. Proper time evolution of cr, N±,T and 

We have determined 5"(m) in terms of the mode functions using two different methods: 
the explicitly renormalized equation, eq. (|77| ) and the combined set of eq. ([7l|) and eq. 



33 




FIG. 9: Evolution o" as a function of u for 3 different values of a. First figure is for & = W ^ , 
second figure is for a = 10~^ and third figure is for a = 10~^. 




0.8 



0.6 



0.6 





FIG. 10: Evolution cr as a function of u for 3 different values of a. First figure is for a = 10 ^ , 
second figure is for a = 10^^ and third figure is for a = 10~^. 
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T and mu on u=1 .0 and 2.5 



2.0 




FIG. 11: Effective temperature and chemical potential as a function of hrj for u=1.0 and u=2.5 

([f6|). The latter set of equations also allows one to check whether one has included enough 
mode functions for the coupling constant to flow logarithmically as in the continuum limit. 
In the previous section we have discussed how these numerical results depended on vari- 
ous discretization parameters as well as the small initial value of the the explicit symmetry 
breaking. From the mode functions, using eq. ( |124| ), one determines the Bogoliubov co- 



efficints and then determines N±{k,u) from eq. (|127|) . By comparing eq. ( |127|) with an 
equilibrium paramatrization: 

N^{k,u) = ll + eM ^^^^^^ ] , (139) 



rp 

where 0;^ = y ^j:^ + we then obtain two equations for the two parameters T{k,u) and 
n{k,u) as a function of A;,,. When these quantities are independent of fc^ this defines the 
proper time evolving temperature and chemical potential. Some indication that an equilib- 
rium parameterization is possible in a 1+1 dimensional field theory evolution in the LOLN 
approximation was already shown in the work of Aarts et. al. |3^. In our simulation we 
find that T and fi are independent of k^^ (except at high momentum) until the system un- 
dergoes the chiral phase transition. A typical example of the dependence of T and fi on /c^ 
is shown in fig. |TI|. For these initial conditions the phase transition occurs near u = 2. The 
connection between kj^ and j of the plot is given in eq. (pH]), namely: k[j] = 27r(j — l/2)/L. 

Because the system goes out of equilibrium after the phase transition, determining T and 
H in that regime is somewhat arbitrary and we use a value averaged over /c^. As we have 
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FIG. 12: Typical evolution of o"/m, fi and T as a function of u for initial conditions below the 
phase transition temperature 

shown earher, if the system evolves in local thermal equilbrium in a massless phase, then 
T falls as 1/r. If there is also chemical equilibrium then fi falls identically to T. We will 
find this is precisely true before the phase transition when there is a second order phase 
transition. For the first order transition fi appears to fall faster than T. We also showed 
that if local chemical and thermal equilbrium are maintained, then the spectrum of particles 
and antiparticles when plotted against the dimensionless momentum /c^ = kr should be 
independent of r. Thus any change in this spectra is an indication of the system going out 
of equilibrium. We expect because of the latent heat released during a first order transition 
that the distortion of the spectra would be greatest in that case which is what we will find 
below. For the time evolution of the single particle distribution, we pick two time u = 1.5 
and u = 2.5 which are before and after the phase transition. 

Let us now focus on the four initial conditions described earlier and look simultaneously 
at the time evolution of the order parameter a, T and /i as a function of u = InnifT. We 
separately plot the evolution of N±{kr^, u). In all these plots all dimensionful parameters are 
scaled by the mass of the Fermion in the broken symmetry vacuum state. 
case (1) fiQ = .2]To = .3 - no phase transition 

Here we start below the phase transition so that the Fermions initially have non-zero mass. 
The results for the order parameter cr as well as for T and /i are displayed in Fig. |l^. We see 
here that although the chemical potential goes to zero as the plasma expands, after u = 1.5, 
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FIG. 13: Evolution of N± as a function of u when there is no phase transiton. The momentum 
displayed is A;^ = kr. 

all the parameters become relatively independent of u and are defined by an equilibrium 
freezeout temperature of approximately .2m f. 

For the broken symmetry case, the plasma essentially stays in equilibrium throughout the 
evolution in that the proper time evolving N± is relatively independent of r and maintains 
its Fermi-Dirac shape. This is seen in Fig. |r^. 
case (2) /io = •5;To = .5 - second order phase transition 

The results for starting in the "unbroken" mode and traversing a second order phase tran- 
sition are shown in Fig. |14[ 

We notice that (j{t) shows a sharp transition during evolution from the unbroken mode 
(a = 0) to the broken symmetry mode. Before the phase transition the temperature falls 
consistent with the equation of state p = e. The chemical potential follows the temperature 
in that regime which means the system is also in chemical equilibrium. After the phase 
transition, the chemical potential goes to zero wheras the temperature freezes out at a 
T ^ .35m f. After the phase transition, there is now a mass scale mf which leads to 
oscillations of a around the final vacuum value. As discussed earlier, to obtain non-trivial 
dynamics we chose a small explicit symmetry breaking parameter cr/mf = 10^"^. 

Going through a second order phase transition does produce a noticable effect in distorting 
the Fermi-Dirac distribution as shown in Fig. ^ 
case (3) Ho = .6; Tq = .32 -traversing the tricritical point 
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FIG. 14: Evolution of T,/x and o" as a function of u for second-order phase transition 
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FIG. 15: Evolution of N± as a function of u when there is a second order phase transiton. The 
momentum displayed is /c^ = kr 



The results for this numerical simulation are shown in Fig. [16. 

We notice that (j{t) again shows a sharp transition during evolution from the unbroken 
mode ((J = 0) to the broken symmetry mode and again the temperature falls consistent with 
the equation of state was p = e before the transition and frezes out again with T ^ .35m f. 
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FIG. 16: Evolution of T,n and a as a function of u passing throught the tricritical point 
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FIG. 17: Evolution of N± as a function of u when one traverses the tricritical point. The momentum 
displayed is krj = kr. 



However the chemical potential in this case falls faster than the temperature. When one 
traverse the tricritical point the distortion of the Fermi-Dirac distribution is greater than 
for the second-order phase transition shown in Fig.|T^. 

case (4) iiQ = .8; To = .3 -first order phase transition 
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FIG. 18: Evolution of T,fx and o" as a function of u passing through first order phase transiton. 
Finally we present results when we start in the unbroken mode and go through a first order 



phase transition. First we display the a, /i, T in Fig. |18[ The results for these variables 
are qualitatively the same as for going through the tricritical point, however the chemical 
potential falls even faster in this case. By comparing the three different cases where there is a 
phase transition, we find among the parameters T,fi and a, only the behavior of the chemical 
potential fi is effected by the order of the phase transition. The chemical potential falls as 
i for the second order phase transition and faster for the first order transition suggesting a 
deviation from local chemical equilibrium. 

The greatest change in the Fermi-Dirac Distribution also occurs for the first order tran- 
sition. This is a result of converting (the small amount of) latent heat into pair production 
and is seen in Fig. W^. 



As we recall from our study of the effective potential, for the GN model the first order 
phase transition is not very strong. The difference in energy density between the false and 
true vacuum (measured in units of rnj'^) is wheras the height of the barrier at phase 
coexistence (seen from Fig. 2) is less than 10 % of this difference. 

B. Numerical results for the energy density and pressure 

In local equilibrium in the massless phase, we have shown that the equation of state 
would he e = p and these quantities would fall as 1/t^ as shown in Fig. ^. In our field 
theory simulations, we will find that this behavior is followed until the phase transition 
occurs. After that the energy density and pressure diverge from each other and oscillate. 
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FIG. 19: Evolution of N± as a function of u when there is a first order phase transiton. 

These oscillations would be damped if we were to go beyond mean field theory and include 
hard scatterings between the Fermions. After the phase transition we find that the energy 
density oscillates around the true broken symmetry values discussed earlier, namely 

1 

'° - "4^ • 

For these simulation we will assume the initial conditions described earlier and plot the 



renormalized energy density and pressure described by eq. (|100|) and eq. (|107|) . Starting 
in the massless (unbroken symmetry) regime we will find p = e and this fall off pertains 
until one goes through the phase transition, after which the system is no longer in local 



equilibrium. In Fig. ED] we show the result of plotting the energy density e{u) and pressure 



p for three initial conditions. 

In Fig. ^ we plot the pressure and energy density for the three initial data (cases 
(2), (3), (4) ) for which there is a phase transition. We see that before the phase transition, p 
tracks e. After the transition, the system goes out of equilibrium and then oscillates about 
its vacuum value because there are no damping mechanisms. 



IX. "PION" CORRELATION FUNCTION 

One interesting question we can ask is how correlation lengths change as we go through 
the phase transition. In our toy model there is no physical pion bound state in our ap- 
proximation, but we still can define an effective field for the pion and study the spatial 
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FIG. 20: Evolution of p and e as a function of u for three initial conditions. 
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FIG. 21: From left to right: evolution of the pressure and energy density as a function of u when 
there is a second order phase transition; evolution of the pressure and energy density as a function 
of u when the trajectory passes through a tricritical point; evolution of the pressure and energy 
density as a function of u when there is a first order phase transition. 

dependence of the pion green's function. Other correlation functions one can study in a 
similar fashion are the density fluctuations as well as the condensate fluctuations. 
We can define an effective (neutral) pion field via 

Ti{x) = cipi{x)i-f^tpi{x) , (140) 

where c is a constant. Using our mode expansion we find that 
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J dk,k,[2 - N+{k) - N.{k)]^\fk\' . (141) 



Because the integrand in eq. (|141|) is odd, the expectation value is zero (otherwise there 



would be spontaneous breakdown of parity). 

For the equal time correlation function in LOLN, we obtain the usual Fermion self energy 
loop. This depends on the time evolving distribution of Fermions and anti-Fermions. Apart 
from an overall constant one can write the connected correlation function in the form: 

D{r] - t(- r) = {'Ki'n, T)n{i, r)), = -jT[j'S{r] - t(- r)j'S{r]' - r^; r)] , (142) 

where at equal times, the propagator is just 

Siv,v';r)af^ = {[Mv,r),4{v',r)]) . (143) 

Here a, [3 take on the values {1,2} and are the spinor indices. Using the mode expansion 
eq. (3.20) we find that the equal time propagator can be written as 

SM;rU = /^e^'''^''-'''^[(l-2iV4A;))0+Jr)0l+(r) 

+ (2iV_(A;)-l)0,Jr)0l^(r)] . (144) 

We could also have use the mode expansion in terms of adiabatic mode functions, eq. (7.1) 
and obtained an expression for the equal time propagator in terms of the time evoloving 
adiabatic number distributions. 
Evaluating the trace we obtain 

r'Div - v'; r)= f ^^e^^'-^^^^-^'^Dik, g; r) , (145) 

where 

D{k, q) = ([1 - 2iV+(A;)][l - 2iV+(g)] + [1 - 2N^ik)][l - 2iV_(g)])Fi(fc, q) 

+ ([1 - 2N+ik)][l - 2iV_(g)] + [1 - 2iV+(g)] .[1 - 2iV_(fc)])F2(fc, q) (146) 

and 

Fi{k, q; r) = \fkf\fgf{{kAg + qA^f + [fc(l^, + a) - qi^k + a)f} , (147) 

and 

F2{k, q; t) = \fk\^\f,\\[kq - A^A, + {n^ + a)(^], + a)]^ 

+ [Afe(f2, + + A^iClk + a)]2) . (148) 
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FIG. 22: Evolution of the Correlation function as a function of time. This is for a first order 
transition. 



If we were solving a 3+1 dimensional 4-fermion model with actual pion composite particles 
we would now be in a position to determine the single particle distribution function for the 
pions from the Wigner distribution functions which is just a particular Fourier transform of 
the Green's function over the relative coordinate. 

In Figs. 0, ^ we plot r^D(?7, Tf) as a function of t] for cases (2) and (4) for a sequence 
of times starting near the onset of the phase transition. 



X. CONCLUSION 



In this paper we have performed simulations in the Gross-Neveu model of an expanding 
plasma of fermions and antifermions. We have chosen inital conditions where the density 
matrix is described by single particle distribution functions pertinent to a plasma initially 
in local thermodynamic and chemical equilibrium. The model was treated in the leading 
order in large N approximation. In this approximation the phase diagram at finite chemical 
potential and temperature shares features with that of massless 2-fiavor QCD. We found 
that if we start in the unbroken symmetry phase, the system remains in equilibrium until the 
phase transition and then goes rapidly out of equilibrium as the phase transition is traversed. 
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FIG. 23: Evolution of the Correlation function as a function of time. This is for a second order 
transition. 



The effects of the phase transition are greatest when we traverse a first order phase transition 
and are most noticable in the anti-fermion distribution function. If these effects survive 
hard scatterings then this should have an effect on the distribution of dileptons, just as the 
overpopulation of soft pions during DCC production effected the distruibution of dileptons as 



discussed by us earlier |^ . We also find that before the phase transition, the system behaves 
identically to an ideal fluid in local thermal equilibrium with equation of state p = e. After 
the phase transition, the system quickly reaches its true broken symmetry vacuum value 
for the fermion mass and for the energy density. Since hard scatterings are ignored in this 
approximation, the competition between the expansion of the plasma and the competing 
process for re-equilbration could not be studied here. Also by restricting our simulations to 
inhomogeneities which are boost invariant, we were not able to look at bubble nucleation. 
In future investigations we will attempt to remedy the shortcomings just mentioned. We 
will perform mean field simulations for inhomogeneous initial conditions which we discussed 
before |^ and which have already been undertaken in scalar field theory in 1+1 dimensions 
|3^ . There are also now two different approaches for going beyond 1/A^ based on Schwinger 
Dyson equations which we hope we can implement to study whether rethermalization can 
occur for the type of expansion expected following a heavy ion collision. These approaches 
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are based on different approximations for the generating functional for the 2-PI irreducible 
graphs ||3^ |^ and should allow us to study the question of rethermalization. We also want 
to extend our simulation to a more realistic 0(4) 4- Fermi Model in 3 + 1 dimensions so that 
we can directly study the effect of the phase transition on pion correlation functions. 
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We will choose for convenience the following representation for the matrices 7*^ and 7 
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7 



Thus the spinor eigenstates of ij^ are 



r 

1 



X 



In terms of this exphcit representation we find that the momentum space wave functions 
have the form: 



<f>k 



(A.2) 
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Correlation function: m = 0.66, T= . 3 at u = 1 
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